AFRL-RV-PS- 

TR-2017-0156 


AFRL-RV-PS- 

TR-2017-0156 


IMPACT  OF  POLARIZING  N ON -LAMBERTIAN 
SURFACE  AND  VOLUME  SCATTERING  ON 
POLARIZED  LIGHT  SIGNATURES:  IMPORTANCE 
TO  REMOTE  SENSING 


Knut  Stamnes  and  Zhenyi  Lin 


Stevens  Institute  of  Technology 
1  Castle  Point  on  Hudson 
Hoboken,  NJ  07030-5906 


08  December  2016 


Final  Report 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  IS  UNLIMITED. 


AIR  FORCE  RESEARCH  LABORATORY 

Space  Vehicles  Directorate 

3550  Aberdeen  Ave  SE 

AIR  FORCE  MATERIEL  COMMAND 

KIRTLAND  AIR  FORCE  BASE,  NM  87117-5776 


DTIC  COPY 


NOTICE  AND  SIGNATURE  PAGE 


Using  Government  drawings,  specifications,  or  other  data  included  in  this  document  for 
any  purpose  other  than  Government  procurement  does  not  in  any  way  obligate  the  U.S. 
Government.  The  fact  that  the  Government  formulated  or  supplied  the  drawings, 
specifications,  or  other  data  does  not  license  the  holder  or  any  other  person  or  corporation; 
or  convey  any  rights  or  permission  to  manufacture,  use,  or  sell  any  patented  invention  that 
may  relate  to  them. 

This  report  is  the  result  of  contracted  fundamental  research  which  is  exempt  from  public 
affairs  security  and  policy  review  in  accordance  with  AFI  61-201,  paragraph  2.3.5. 1.  This 
report  is  available  to  the  general  public,  including  foreign  nationals.  Copies  may  be  obtained 
from  the  Defense  Technical  Information  Center  (DTIC)  (http://www.dtic.mil). 


AFRL-RV-PS-TR-20 17-0 156  HAS  BEEN  REVIEWED  AND  IS  APPROVED  FOR 
PUBLICATION  IN  ACCORDANCE  WITH  ASSIGNED  DISTRIBUTION  STATEMENT. 


//SIGNED// 


//SIGNED// 


Jeannette  van  den  Bosch 
Program  Manager,  AFRL/RVBYI 


Dr.  Thomas  R.  Caudill,  Acting  Chief 
AFRL  Battlespace  Environment  Division 


This  report  is  published  in  the  interest  of  scientific  and  technical  information  exchange,  and  its 
publication  does  not  constitute  the  Government’s  approval  or  disapproval  of  its  ideas  or  findings. 


Approved  for  public  release;  distribution  is  unlimited. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE  3.  DATES  COVERED  (From  -  To) 

08-12-2016  Final  Report  09  Jun  2015  -  01  Oct  2016 


4.  TITLE  AND  SUBTITLE 

Impact  of  Polarizing  Non-Lambertian  Surface  and  Volume  Scattering  on  Polarized  Light 
Signatures:  Importance  to  Remote  Sensing 


6.  AUTHOR(S) 

Knut  Stamnes  and  Zhenyi  Lin 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 

FA9453- 15- 1-0076 


5c.  PROGRAM  ELEMENT  NUMBER 

62601F 


5d.  PROJECT  NUMBER 

1010 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Stevens  Institute  of  Technology 
1  Castle  Point  on  Hudson 
Hoboken,  NJ  07030-5906 


5e.  TASK  NUMBER 

PPM00026607 


5f.  WORK  UNIT  NUMBER 

EF129628 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory 
Space  Vehicles  Directorate 
3550  Aberdeen  Avenue  SE 
Kirtland  AFB,  NM  87117-5776 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFRL/RVBYI 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

AFRL-R  V -PS  -TR-20 17-0156 


14.  ABSTRACT 

A  modernized  and  upgraded  version  of  the  polarized  (vector)  radiative  transfer  algorithm  VDISORT-OD5  will  be  used  to  study  the 
interaction  of  polarized  light  with  realistic  (rough)  surfaces.  The  end  goal  is  to  analyze  the  impact  of  rough  surfaces  on  polarized  light  at 
arbitrary  polar  angles,  and  to  carry  out  a  sensitivity  analysis  of  the  influence  of  rough  surfaces  on  polarized  light.  Although  thermal 
radiation  is  initially  unpolarized  when  emitted  in  the  atmosphere,  it  becomes  polarized  by  scattering  processes  with  aerosols  and  clouds  and 
through  interactions  with  polarizing  surfaces.  The  proposed  work  will  analyze  the  impact  of  rough  surface  and  volume  scattering  on 
polarized  light  at  arbitrary  polar  angles,  and  perform  a  sensitivity  analysis  of  the  influence  of  rough  surface  and  volume  scattering  on 
polarized  light.  Utilization  of  the  polarized  bi-directional  reflection  capabilities  in  VDISORT-MOD5  will  enhance  our  understanding  of  the 
interaction  of  polarized  light  with  realistic  surfaces  and  is  of  utmost  importance  in  remote  sensing  analysis  for  correct  characterization  of 
rough  surfaces  with  respect  to  fundamental  spectral  signature  development,  surface  modeling,  experimental  validation,  synthetic  image 
generation,  and  improvement  of  both  quantitative  remote  sensing  algorithm  performance  and  accuracy  assessment  of  image  analysis. 


15.  SUBJECT  TERMS 

polarized  (vector)  radiative  transfer,  phase  matrix,  phase  function,  Stokes  parameters,  non-Lambertian  reflecting  surfaces 


16.  SECURITY  CLASSIFICATION  OF: 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

c.  THIS  PAGE 

Unclassified 

17.  LIMITATION 

18.  NUMBER 

OF  ABSTRACT 

OF  PAGES 

Unlimited 

40 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Jeannette  van  den  Bosch 


19b.  TELEPHONE  NUMBER  (include  area 
code) 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  239.18 


This  page  is  intentionally  left  blank. 


Approved  for  public  release;  distribution  is  unlimited. 


Table  of  Contents 


1.  Summary . 1 

2.  Upgrade  of  VDISORT . 1 

2.1  General  discussion  of  arbitrary  polar  user-angle  output . 1 

2.2  Analytic  representation  of  the  source  function . 3 

2.3  Upgrade  of  the  arbitrary  angle  output . 4 

2.4  Validation  of  the  arbitrary  angle  output  -  black  surface . 4 

2.5  Validation  of  user-angle  output  -  wind-roughened  ocean  surface . 6 

3.  Sensitivity  Study . 6 

4.  Description  of  the  code . 11 

4.1  Calling  VDISORT:  benchmark  computations  and  other  tests . 11 

4.2  Model  setup . 11 

5.  References . 12 

6.  Appendix:  Theoretical  Background . 13 

6.1  RTE  for  polarized  radiation . 14 

6.2  Solution  of  the  RTE . 16 

6.3  Discrete  ordinate  equations  -  Compact  matrix  formulation . 21 

6.4  Discrete  ordinate  solutions . 24 


Approved  for  public  release;  distribution  is  unlimited. 


1 


List  of  Figures 


Figure  1:  Comparison  of  output  for  (i)  quadrature  angles;  (ii)  spline  interpolation;  and  (iii) 
arbitrary  angle  output  ('analytic') . 2 

Figure  2:  Kokhanovsky's  benchmark,  I  component . 4 

Figure  3:  Kokhanovsky's  benchmark,  Q  component . 5 

Figure  4:  Kokhanovsky's  benchmark,  U  component  (left)  and  V  component  (right) . 5 

Figure  5:  Comparison  of  SeaDAS  Rayleigh  tables,  VDISORT  Rayleigh  tables  (vector),  DIS¬ 
ORT  Rayleigh  tables  (scalar),  and  VDISORT  full  radiances  (vector)  for  a  relative  azimuth 
angle  of  0° . 7 

Figure  6:  Comparison  of  SeaDAS  Rayleigh  tables,  VDISORT  Rayleigh  tables  (vector),  DIS¬ 
ORT  Rayleigh  tables  (scalar),  and  VDISORT  full  radiances  (vector)  for  a  relative  azimuth 
angle  of  90° . 7 

Figure  7:  Comparison  of  SeaDAS  Rayleigh  tables,  VDISORT  Rayleigh  tables  (vector),  DIS¬ 
ORT  Rayleigh  tables  (scalar),  and  VDISORT  full  radiances  (vector)  for  a  relative  azimuth 
angle  of  180° . 8 

Figure  8:  Test  of  Stream  number:  I  component . 8 

Figure  9:  Test  of  Stream  number:  Q  component . 9 

Figure  10:  Test  of  Stream  number:  U  component  (left)  and  V  component  (right) . 9 

Figure  1 1 :  Comparison  of  I  component  w/  and  w/o  DeltaFit  for  16  streams . 10 

Figure  12:  Comparison  of  I  component  w/  and  w/o  DeltaFit  for  32  streams . 10 


Approved  for  public  release;  distribution  is  unlimited. 


11 


Final  Performance  Report 

Award  FA9453-15-1-0076 

Impact  of  polarizing  non-Lambertian  surface 
and  volume  scattering  on  polarized  light 
signatures:  Importance  to  remote  sensing 

Knut  Stamnes,  PI 
Stevens  Institute  of  Technology 


1.  Summary 

In  2015  -  2016,  we  have  focused  on  upgrading  VDISORT  to  improve  the 
computation  of  the  Stokes  components  at  arbitrary  polar  viewing  angles.  We 
also  have  carried  out  sensitivity  tests  for  non-Lambertian  surfaces.  The  work 
carried  out  includes  the  following  accomplishments: 

•  An  upgrade  of  VDISORT  to  allow  for  output  at  arbitrary  (user)  polar 
viewing  angles  for  the  3- vector  (I,  Q ,  U)T  as  well  as  the  4- vector  (I,  Q,  U,  V)T 
Stokes  vector  has  been  completed.  Here  T  denotes  the  transpose. 

•  An  extension  from  a  single  layer  to  a  general  multiple-layer  system  has 
been  completed  for  the  user  polar  viewing  angle  output. 

•  The  user  polar  viewing  angles  output  has  been  extended  from  working 
for  a  black  surface  only  to  a  partially  reflecting  surface  represented  by 
Lambertian  as  well  as  general  non-Lambertian  surfaces. 

•  Sensitivity  studies  of  how  the  user  polar  angle  output  depends  on  the 
number  of  streams  have  been  carried  out. 


2.  Upgrade  of  VDISORT 

2.1.  General  discussion  of  arbitrary  polar  user-angle  output 

The  discrete  ordinate  method  provides  Stokes  vector  components  at  discrete 
polar  angles,  called  the  quadrature  angles.  However,  many  practical  problems 
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including  remote  sensing  applications  require  Stokes  vector  components  at  sen¬ 
sor  observing  angles.  Even  though  one  could  envision  generating  results  at 
sensor  observing  angles  by  interpolating  or  extrapolating  the  quadrature  values 
to  such  angles,  accurate  interpolation  of  these  values  is  difficult,  particularly  for 
optical  depths  close  to  the  top  boundary.  For  example,  Fig.  1  shows  an  example 
of  inaccurate  interpolation  (’spline’)  for  a  benchmark  case 
1989)  that  has  been  reproduced  by  VDISORT.  We  note  that  the  ‘spline’  interpo- 


Figure  1:  Comparison  of  output  for  (i)  quadrature  angles;  (ii)  spline  interpolation;  and  (iii) 
arbitrary  angle  output  (‘analytic’). 

lation  produces  large  oscillations,  wheres  the  integration  of  the  source-function 
technique,  see  Eqs.  (1)  and  (2),  give  accurate  (‘analytic’)  results  that  agree  with 
the  benchmark  values. 

The  upgrade  of  VDISORT  to  give  output  at  user-specified  polar  angles,  pro¬ 
vides  a  robust  and  accurate  way  to  generate  results  at  arbitrary  polar  angles 
analytically.  For  a  plane  parallel  atmosphere,  the  Stokes  components  at  an  arbi¬ 
trary  polar  angle,  9  =  cos-1  /x,  can  be  computed  by  an  integration  of  the  source 
function.  In  plane-parallel  geometry  the  Stokes  vector  I(r,  /x,  </>)  =  (I,  Q,  U,  V)T, 
where  T  denotes  the  transpose,  can  be  expanded  in  a  Fourier  series  in  a  Fourier 
series  (see  Eq.  (22)  of  the  Appendix),  and  the  Fourier  components  can  be  ex¬ 
pressed  as 


C  +(r,9) 

=  C+(r6,/x)e-^-^^  + 

r  -s™+(t,/x)e-(t-T)/'i 

(1) 

lr  M 

I 

=  ir(T0,M)e-(T-To)Ai  + 

(2) 

>T0  H 

at  an  arbitrary  angle  0  =  cos _1/i  in  the  downward  [I™“(r,  /x)]  and  upward 
[I ™+(r, /x))]  directions.  Here  r0  =  0  is  the  optical  depth  at  the  top  of  the  slab, 
and  77,  is  the  optical  thickness  of  the  slab. 

A  very  important  reason  to  develop  an  analytic  arbitrary  polar  angle  output 
capability  in  VDISORT  is  that  MODTRAN  makes  use  of  this  feature,  which  is 
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particularly  useful  when  considering  spherical  refractive  geometry.  Computation 
of  Im±(r,  p)  at  arbitrary  values  of  cosine  of  the  polar  angle  p  based  on  Eqs.  (1) 
and  (2)  is  similar  to  that  used  in  the  scalar  radiative  transfer  model  (DISORT) 
which  has  already  been  implemented  into  MODTRAN  (see  MODTRAN  ATBD 
for  more  details).  However,  the  implementation  of  VDISORT  into  MODTRAN 
becomes  more  complicated  because  (i)  the  four  Stokes  components  are  coupled 
and  must  be  computed  together;  (ii)  the  introduction  of  the  ‘V’  component 
produces  some  complex  eigenvalues  and  eigenvectors  which  require  additional 
attention  and  processing.  To  be  more  specific,  we  may  define  the  incremental 
Stokes  vector  for  a  line-of-sight  (LOS)  segment  using  Eqs.  (1)  and  (2): 

rTh 

A  C+(r,/z)  =  C+(T,/z)-C+(Tb,M)e-^-^=  /  -S^+(t,  p)e~^^ 

J  T  lL 

(3) 

CT  (ii 

Air(r,M)  =  ir(T,AO-ir(^o,M)e-(T-To)/^=  /  -S™-(t,/1)e-{T-‘)/'‘. 

J  To  A* 

(4) 

For  a  single-layer,  homogeneous  slab,  the  LOS  segment  will  be  a  straight  line, 
while  for  a  multi-layer,  plane- parallel  atmosphere,  each  layer  LOS  corresponds  to 
a  straight  line  segment.  As  explained  in  some  detail  in  the  MODTRAN6  ATBD, 
in  a  realistic  curved  (and  refractive)  atmosphere,  each  layer  should  be  divided 
into  several  segments,  and  the  value  of  p,  will  be  different  in  each  segment  along 
the  curved  LOS.  As  in  the  scalar  DISORT  case,  implementation  of  arbitrary 
user-angle  output  greatly  facilitates  the  computation  of  LOS  segment  Stokes 
components,  and  the  incorporation  of  VDISORT  into  MODTRAN  in  spherical, 
refractive  geometry. 

2.2.  Analytic  representation  of  the  source  function 

As  shown  in  Section  6.1,  the  source  function  can  be  written  as 

Sq  (t,u)  =  S™(r,u)  +  Q™(t,m)  a  =  cors  (5) 

where  [see  Eq.  (28)] 

Q  =  X.™(T,u)e~T^°+  6ca50m[l-w(T)]St(T)  (6) 

and  the  contributions  due  to  multiple  scattering  are  given  by 

j  =  —  N  { 

3  t^O 

-  P r(r,  Uj,u)  I r(r,  «,■)}  -  Q T(r,  u)  (7) 

s ?(T,u)  =  ^  E  T(r,Uj) 

j  =  —  N  y. 

3  t^O 

+  P T(r,  ujtu)  I?(t,  Uj-)}  -  Q?{t,  u).  (8) 
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We  see  that  the  source  function  is  obtained  from  a  summation  of  quadrature 
angle  output  where  fij  indicates  the  jth  of  the  2 N  quadrature  angles.  For 
instance,  I™(t,  /Uj, )  is  the  Stokes  vector  evaluated  at  the  jth  quadrature  angle 

li  j. 

2.3.  Upgrade  of  the  arbitrary  angle  output 

The  upgrade  of  arbitrary  angle  output  is  based  on  the  progress  in  the  pre¬ 
vious  year.  In  2015,  the  arbitrary  angle  output  could  only  be  computed  under 
a  black  lower  boundary,  a  single  layer  and  for  the  first  three  ([/,  Q ,  U]T)  Stokes 
components.  Now  in  2016,  the  arbitrary  angle  output  in  VDISORT  is  improved 
to  provide  the  support  for:  (i)  the  computation  of  the  full  Stokes  vector  output, 
including  the  lV’  component;  (ii)  a  multi-layer  system;  (iii)  a  reflectance  from 
a  general  (Non-Lambertian)  lower  boundary. 

2.f.  Validation  of  the  arbitrary  angle  output  -  black  surface 

To  validate  the  arbitrary  user-angle  output,  we  compared  with  Kokhanovsky’s 
benchmark  computations  of  transmittances  and  re¬ 
flectances  of  Stokes  vector  components  for  a  single  medium  consisting  of  ei¬ 
ther  molecules  only,  aerosols  only,  or  cloud  droplets  only.  The  aerosol  and 
cloud  particles  were  considered  to  be  spheres  whose  Inherent  Optical  Proper¬ 
ties  (IOPs)  were  computed  by  Mie-theory.  In  the  previous  year,  we  had  tested 
Kokhanovsky’s  benchmark  in  the  quadrature  angles  and  obtained  a  good  agree¬ 
ment.  This  year  we  tested  the  arbitrary  angle  output. 


Figure  2:  Kokhanovsky’s  benchmark,  I  component. 


Figures  2-4  show  the  comparison  of  benchmark  computations  for  an  aerosol 
case  with  VDISORT  results  for  200  polar  output  angles.  In  this  test  we  used  192 
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O  vdisort  user  angle 
- benchmark-aerosol-refl 

/ 

Figure  3:  Kokhanovsky’s  benchmark,  Q  component. 


Figure  4:  Kokhanovsky’s  benchmark,  U  component  (left)  and  V  component  (right). 


streams  and  Delta-M  scaling  (Wiscombc.  1971  .  The  200  user-defined  angles 
are  shown  as  blue  circles  in  each  figure  with  cos  8  varying  from  —1.0  to  1.0  in 
steps  of  0.01.  In  order  to  test  of  the  user-angle  output  for  a  multi-layer  system, 
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we  split  the  single  layer  into  5  identical  layers  with  the  same  optical  properties 
and  the  same  total  optical  thickness  as  for  the  single  layer.  Again,  very  good 
agreement  was  obtained  for  the  user-angle  output. 


2.5.  Validation  of  user- angle  output  -  wind-roughened  ocean  surface 

A  rough  ocean  surface  has  also  been  added  to  VDISORT  and  used  to  test 
the  user-angle  output  implementation.  An  explicit  expression  for  the  reflectance 
matrix  Rrs  at  the  rough  surface  interface  is  given  by 


i  IM4 


Rrs(^S3  4* si  4i) 


cos9s  4|fcj  x  ks\4kfdz 


The  factor 


is  the  probability  density  of  surface  slopes  resulting  from  the  assumption  of  a 
Gaussian  rough  surface ,  and  a2  is  the  mean  square  surface  slope. 

To  test  the  implementation  of  the  rough  surface  in  VDISORT,  we  decided 
to  reproduce  the  SeaDAS  Rayleigh  tables  .In  these 

tables  the  TOA  radiance  contributed  from  a  clear  atmosphere  (molecular  scat¬ 
tering  only)  overlying  a  wind-roughened  water  surface  are  tabulated.  In  VDIS¬ 
ORT  we  used  a  Rayleigh  scattering  atmospheric  layer  above  a  rough  surface 
to  reproduce  the  radiances  provided  in  these  SeaDAS  tables  for  a  wind  speed 
of  4.21  m/s.  We  turned  off  the  direct  beam  reflectance,  but  kept  the  diffuse 
reflectance  in  order  to  remove  the  radiance  due  to  sunglint. 

Figures  5-7  show  comparisons  of  radiances  at  412  nm  for  three  different  az¬ 
imuth  angles,  0°,90°,  and  180°,  and  for  a  20°  solar  zenith  angle.  In  general, 
there  is  a  good  match  between  VDISORT  results  and  the  SeaDAS  tables  with 
a  less  than  3%  difference.  By  contrast,  a  difference  up  to  10%  is  found  between 
the  result  produced  by  DISORT  (green  line)  and  the  SeaDAS  tables,  demon¬ 
strating  the  importance  of  polarization  effects.  The  full  radiance  (the  blue  line), 
including  both  the  diffuse  (skyglint)  and  the  direct  beam  (sunglint)  reflectance, 
is  shown  to  illustrate  how  strong  the  sunglint  is  at  different  azimuth  angles. 

3.  Sensitivity  Study 

The  phase  matrix  for  spherical  particles  is  expanded  into  four  independent 
Greek  constants  a\  =  0%  a 3  =  oq,  /3\,  fa-  In  Kokhanovsky’s  benchmark  for 
aerosol  particles,  each  constant  is  represented  by  more  than  900  terms,  which 
implies  that  900  streams  would  be  needed  to  obtain  very  accurate  results.  In 
the  last  section,  we  used  192  streams  as  well  as  the  Delta-M  method 
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0.1 


channel  =  412nm,  SZA  =20  °,  wind  speed  =  4.2106  m/s,  A</>  =  0° 


Senor  Zenith  Angle 


Figure  5:  Comparison  of  SeaDAS  Rayleigh  tables,  VDISORT  Rayleigh  tables  (vector),  DIS¬ 
ORT  Rayleigh  tables  (scalar),  and  VDISORT  full  radiances  (vector)  for  a  relative  azimuth 
angle  of  0°. 
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Figure  6:  Comparison  of  SeaDAS  Rayleigh  tables,  VDISORT  Rayleigh  tables  (vector),  DIS¬ 
ORT  Rayleigh  tables  (scalar),  and  VDISORT  full  radiances  (vector)  for  a  relative  azimuth 
angle  of  90°. 


1977  to  match  the  benchmark  results  and  obtained  good  agreement.  However, 
in  practice  we  generally  want  to  use  much  fewer  streams  (e.g.  less  than  32 
streams)  because  the  computational  expense  increases  sharply  as  the  number 
of  streams  is  increased.  To  examine  how  the  accuracy  depends  on  the  number 
of  streams,  we  re-computed  the  aerosol  benchmark  case  using  8,  16,  and  32 
streams,  and  compared  the  results  with  the  accurate  19-stream  results.  These 
comparisons  are  shown  in  Fig.  8  -  10.  We  note  that  the  8-  and  16-stream  cases 
(green  and  blue  curves)  failed  to  match  the  192-stream  results  (red  curves).  The 
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channel  =  412nm,  SZA  =20°,  wind  speed  =  4.2106  m/s,  A <j>  -  180 


Senor  Zenith  Angle 


Figure  7:  Comparison  of  SeaDAS  Rayleigh  tables,  VDISORT  Rayleigh  tables  (vector),  DIS¬ 
ORT  Rayleigh  tables  (scalar),  and  VDISORT  full  radiances  (vector)  for  a  relative  azimuth 
angle  of  180°. 


Figure  8:  Test  of  Stream  number:  I  component. 


32-stream  results  (black  curves)  are  close  to  the  192-stream  results,  except  for 
the  discrepancies  around  the  forward  scattering  (left  bottom  subplot,  around 
60°)  and  backward  scattering  angles  (right  top  subplot,  around  60°)  as  well  as 
some  oscillations. 

To  further  improve  the  32-stream  results,  we  introduced  the  Delta-Fit  method 
(Hu  et  al.,  2000),  which  provides  a  stronger  truncation  of  the  forward  scattering 
peak,  compared  with  Delta-M  method.  Specifically,  Delta-Fit  method  applies 
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Figure  9:  Test  of  Stream  number:  Q  component. 


Figure  10:  Test  of  Stream  number:  U  component  (left)  and  V  component  (right). 


a  Singular  Value  Decomposition  (SVD)  followed  by  a  Least  Squares  Fitting  to 
the  phase  matrix  to  improve  the  performance  for  small  stream  numbers. 

To  examine  the  effect  of  using  Delta-Fit,  Figures  11  and  12  show  the  im¬ 
provement  of  the  /  component  (radiance)  resulting  from  the  use  of  the  Delta-Fit 
method  for  16  and  32  streams.  The  results  of  Delta-Fit  (red  curves)  show  an 
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improvement  and  better  match  with  the  192  stream  case  (blue  curves).  The 
oscillations  in  32-stream  results  without  Delta-Fit  scaling  (green  curves)  have 
also  been  removed. 


View  Zenith  Angle 


Figure  11:  Comparison  of  I  component  w/  and  w/o  Delt aF it  for  16  streams. 


View  Zenith  Angle 


Figure  12:  Comparison  of  I  component  w/  and  w/o  DeltaFit  for  32  streams. 
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4.  Description  of  the  code 


In  this  section,  examples  are  shown  of  how  to  call  VDISORT  for  different 
conditions. 

4-.1.  Calling  VDISORT:  benchmark  computations  and  other  tests 

Table  shows  a  list  of  benchmark  tests  for  VDISORT.  These  benchmarks  are 
from  Garcia  and  Siewert  (1989),  Siewert  (2000),  and  Kokhanovsky  et  al.  (2010). 


Script  file 

Description 

Matlab  plots 

r  un_kokh2  0 1 0  _r  ay  leigh .  sh 
r  un_kokh2  010  -aerosol .  sh 
r  un_kokh2  0 1 0  _cloud .  sh 
r  un_siewert  1 9 89 .  sh 
run_siewert2000 .  sh 

Kokhanovsky’s  Rayleigh  benchmark 
Kokhanovsky’s  Aerosol  benchmark 
Kokhanovsky’s  Cloud  benchmark 
Garcia  and  Siewert  (1989)  benchmark 
Siewert  (2000)  benchmark 

benchmark_vdisort_kokh_rayleigh.m 
benchmark_vdisort_kokh_aerosol.m 
benchmark_vdisort_kokh_cloud.m 
benchmark_siewert 1 989 .  m 
benchmark_siewert2000.m 

To  run  the  benchmarks,  simply  execute  the  corresponding  scripts  in  the  table. 
To  verify  that  the  results  are  correct,  a  Matlab  code,  located  in  the  directory 
‘./plots’,  is  provided  to  produce  plots  that  compare  VDISORT  output  with  the 
benchmark  data. 

Similarly,  there  are  some  other  tests  for  the  arbitrary  angle  output,  includ¬ 
ing  multiple  (2)  layers  as  well  as  for  a  general  Non-Lambertian  surface  shown 
in  table. 


Script  file 

Description 

run_usrang_21ayers-kokh2010_aerosol.sh 
run_usr  ang_21ayer  s_siewert  1 989 .  sh 
run_usrang_21ayers_siewert2000.sh 
run_usrang_vacuum_cmbrdf 
run usr  ang r  ayleigh-cmbrdf 

Arbitrary  angle  and  split  2  layers  for  Kokhanovsky’s  Aerosol  benchmark 
Arbitrary  angle  and  split  2  layers  for  Garcia  and  Siewert ’s  benchmark 
Arbitrary  angle  and  split  2  layers  for  Siewert  benchmark  in  2000 

Arbitrary  angle  for  the  Cox-Munk  rough  surface 

Arbitrary  angle  for  a  Rayleigh  layer  on  top  of  the  Cox-Munk  rough  surface 

Finally,  there  are  script  tests  included  in  table  that  can  be  used  to  reproduce 
some  DISORT  3  tests.  To  reproduce  the  results  from  the  scalar  model,  all  ele¬ 
ments  of  the  four  by  four  scattering  phase  matrix,  except  the  first  one,  are  set 
to  zero. 


Script  file 

Description 

run_disort3_case2a.sh 
run_disort3_case2b.sh 
run_disort3_case2c.sh 
run_disort3_case2d.sh 
run_disort3_case6c.sh 
r  un_disort  3_case  1 4b .  sh 

Reproduce  DISORT  test  case  2a:  Rayleigh  scattering,  Beam  source 
Reproduce  DISORT  test  case  2b:  Rayleigh  scattering,  Beam  source 
Reproduce  DISORT  test  case  2c:  Rayleigh  scattering,  Beam  source 
Reproduce  DISORT  test  case  2d:  Rayleigh  scattering,  Beam  source 
Reproduce  DISORT  test  case  6c:  Lambert  surface 

Reproduce  DISORT  test  case  14b:  Rough  sea  surface  (Non-Lambertian) 

4-2.  Model  setup 

To  build  the  personal  settings,  one  may  need  to  know  model  setups.  Below 
shows  some  of  important  Model  inputs: 
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Name 

Description 

Lambert 

True  =  lambert  surface;  False  =  general  surface 

ALBEDO 

if  Lambert  =  true,  it  gives  the  value  of  albedo 

NLYR 

number  of  layers 

NLEG 

number  of  GSF  terms  that  expands  phase  matrix 

NAZ 

NAZ  =  NLEG  -1 

NSTR 

NSTR  =  NLEG 

USRANG 

True  =  output  in  arbitrary  angles;  False  =  output  in  quadrature  angles 

NUMU 

if  USRANG  =  true,  it  gives  the  number  of  arbitrary  polar  angles 

UMU 

if  USRANG  =  true,  it  gives  the  cosine  of  arbitrary  polar  angles  (1:NUMU) 

USRTAU 

True  =  specify  user  output  optical  thickness;  False  =  layer  optical  thickness 

NUTAU 

if  USRTAU  =  true,  it  gives  the  number  of  output  optical  thickness 

UTAU 

if  USRTAU  =  true,  it  gives  the  optical  thickness  for  output  (1:NUTAU) 

NUPHI 

the  number  of  output  azimuth  angles 

UPHI 

azimuth  angles  for  output  (1:NUPHI) 

UMUO 

cosine  of  Solar  Zenith  Angle 

PHIO 

solar  azimuth  angle  (PHIO  =  0) 

FBEAM 

FBEAM  =  so  (solar  incidence) 

And  the  following  table  shows  some  of  important  model  parameters,  which  are 
located  in  the  directory  ‘/src/PARM_F90.h’.  It  is  worth  mentioning  that  once 
changes  have  been  made  to  the  input  constants,  it  is  then  necessary  to  clean 
all  binary  files  and  re-compile  the  whole  model.  This  “clean-up”  can  be  done 
by  simply  going  to  the  “/src”  directory  and  typing  the  command  “gmake  clean” . 


Variable  name 

Description 

Reduce 

Reduce=‘l’:  apply  reduction  of  order;  (fast) 

Reduce=‘0’:  does  not  apply  reduction  of  order  (slow) 

Recommend  to  set  Reduce=‘l’  to  speed  up 

Complx 

Complx=T’:  compute  complex  solution  [/,  Q ,  U,  V]1  ’; 

Complx=‘0’:  ignore  V  component  and  compute  solution  [J,  Q,  U,  0]T 
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6.  Appendix:  Theoretical  Background 

The  radiative  transfer  equation  (RTE)  for  the  total  unpolarized  radiance, 
which  includes  both  thermal  emission  and  multiple  scattering,  can  be  written 

dIu(T,  tt)  _  _j^T,  O)  4.  [1  _  m(v)\Bu(T)  +  fi')  (10) 

drs  47t  _/47r 
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where  Iv(t,  fi)  is  the  radiance,  cu(z'),  the  single-scattering  albedo,  BU(T)  the 
Planck  function,  and  p(Of,  fl)  the  scattering  phase  function.  For  a  plane-parallel 
medium  in  which  the  IOPs  are  assumed  to  vary  only  in  the  vertical  direction 
denoted  by  z,  increasing  upwards,  so  that  the  corresponding  vertical  optical 
depth,  denoted  by  r(z),  is  defined  by 

/»oo 

t(z)=  [a(z')  +  /3{z')]dz'  (11) 

J  Z 

and  hence 

dr(z)  =  —  [a(r)  +  P{r)]dz  (12) 

where  the  minus  sign  indicates  that  r  increases  in  the  downward  direction, 
whereas  z  increases  in  the  upward  direction.  Equation  (10)  pertains  to  the  total 
radiation  field.  In  a  plane-parallel  geometry,  with  a  collimated  beam  incident  at 
the  top  of  the  slab,  one  may  invoke  the  usual  diffuse-direct  splitting  to  obtain 
for  the  diffuse  component  of  the  radiation  field  (omitting  the  subscript  v  and 
the  functional  dependence  on  v  for  notational  convenience) 

u  ^  =  J(r>  u’<t>)  ~  S(T>  u,  <t>)  (13) 

where 


S{t,u,4>)  =  Sb(T,u,<j))  +  [1  -  w(t)\B(t) 

2tt  1 

+  ~^  J  J  p(T,u',<j)'-,u,<j))I(T,u',<j)')du.  (14) 

0  -1 

Here  u  is  the  cosine  of  the  polar  angle  9 ,  (j)  is  the  azimuth  angle,  vj(t)  = 

/3(r)/[a(r)  +  /3(r)]  is  the  single-scattering  albedo,  p(r,  u' ,  <j)’\ u,  <f>)  is  the  scatter¬ 
ing  phase  function,  and  B(t)  is  the  thermal  radiation  field  given  by  the  Planck 
function.  The  term  Sb(r,  u,  <j>)  in  Eci.  (14)  is  the  single-scattering  source  function 
due  to  a  collimated  beam  incident  at  the  top  of  the  upper  slab,  to  be  specified 
below. 

6.1.  RTE  for  polarized  radiation 

To  generalize  Eq.  (13)  to  apply  to  polarized  radiation,  we  note  that  the 
multiple  scattering  term  in  Eq.  (14)  must  be  replaced  by 

2tt  1 

S(r,  u,  (j))  =  J  d4>  J  P(t,  u',  <f>\  u,  cf>)l (r,  v! ,  <j>')du  (15) 

0  -1 

where  I(r,  u' ,  (j)’)  is  the  Stokes  vector,  and  P(r,  u1,  u,  (j))  is  the  scattering  phase 
matrix.  The  first  element  of  the  vector  S(r,  u,  <j>)  represents  the  energy  per  unit 
solid  angle,  per  unit  frequency  interval,  and  per  unit  time  that  is  scattered  by  a 
unit  volume  in  the  direction  ( 9,(f> ).  Hence,  in  a  plane-parallel  (slab)  geometry, 
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the  integro-differential  equation  for  polarized  radiative  transfer  is  expressed  in 
terms  of  a  Stokes  vector  I(r,  u,  <j>)  as 


dl(r,  u,  4>) 

U  J 

where  the  source  function  is 


I(r,u,</>)  -  S(t,u,  4>) 


(16) 


S(t,  u,  <f>) 


j  ^  [  d(j)'  [  dv!  P(r,  u' ,  <f>\ u,  u1,  <j>') 

47r  Jo  J-i 

+  Q  {t,u,4>).  (17) 


The  term  Q(r,  u,  (f>),  due  to  thermal  and  beam  sources  is  given  by: 

Q{t,u,  <j>)  —  P(r’  —^0;  '/’o;  u,  4>)Sbe~T^°  +  [1  -  w(t)\  St(r).  (18) 

The  first  term  on  the  right  hand  side  of  Eq.  (18)  describes  the  incident  beam 
Sf ,  in  direction  (—poJo),  which  is  attenuated  at  depth  r  by  a  factor  e~T^° 
and  undergoes  single  scattering  into  the  direction  (it,  f>) .  For  an  unpolarized 
incident  beam  S/,  has  the  form 

Sb  =  [Jo/2,  Jo/2,  0,  0]T  or  [J0, 0, 0,  0]T  (19) 


where  the  first  or  second  expression  corresponds  to  the  choice  of  Stokes  vector 
representation,  [I\\,  I±,U,V]T  or  [I,Q,U,  V]T .  The  second  term  on  the  right 
hand  side  of  Ecp  (18)  is  due  to  thermal  emission,  which  is  unpolarized  and  given 

by 

St(r)  =  [B(T(T))/2,B(T(r))/2,0,0]T  or  [B(T(r)),  0,  0,  0]T  (20) 

where  B  is  the  Planck  function,  and  where  the  first  or  second  expression  corre¬ 
sponds  to  the  choice  of  Stokes  vector  representation.  We  have  set  p o  =  |uo|  = 
|  cos  1 1  where  9q  is  the  polar  angle  of  the  incident  light  beam. 


6.1.1.  Isolation  of  azimuth  dependence 

We  start  by  expanding  the  phase  matrix  in  a  Fourier  series: 

M 

P(it',  u;  (/>' —  <f>)  =  ^  {P™ (u',u)  cos  m(<j)'  —  (f)+  P™  (t/,  it)  sin  rn(^/  —  (/>)}.  (21) 

m= 0 

To  isolate  the  azimuth  dependence  of  the  radiation  field  we  expand  the  Stokes 
vector  I(r,  u,  4>)  in  Eq.  (16)  and  the  source  term  Q(r,  u,  (j> )  in  Eq.  (18)  in  a  Fourier 
series  in  a  manner  similar  to  the  expansion  of  the  phase  matrix  in  Eq.  (21): 

M 

I  {t,u,4>)  =  Y 

m= 0 
M 

Q  {t,u,4>)  =  Y 

m— 0 


Q™  (r,  u)  cos  m((j) 0  —  </>)+  Q™  (t>  u )  sin  —  </>)  f  (23) 


I ™{t,u)  cos  m(<j)Q  —  </>)+  I™(t,  u)  sin  m(4>0  -  ft) 


(22) 
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where  the  subscript  s  or  c  denotes  sine  or  cosine  mode.  Substitution  of  Eqs.  (21) 


into  Eq.  (18)  and  comparison  with  Eq.  (23)  yields: 

Q  =  X™(T,u)e~r^°+6om[l-w(T)]St(T)  (24) 

CTO-.u)  =  X™(r,W)e-^°  (25) 

where 

X™(r,u)  =  (26) 

XT(t,u)  =  P™(T) — MOi  u)  Sb  (27) 

which  can  be  written  more  compactly  as  (letting  a  =  c,  s  represent  the  cosine 
or  sine  mode): 

Qa(r,u)  =  X™(t,  u)  e~T^°  +  5ca  <50m  [1  —  w(r )]  St(r )  (28) 

where  Sca  =  1  if  a  =  c,  and  Sca  =  0  if  a  ^  c. 


By  substituting  Eqs.  (21),  (22),  and  (23)  into  Eqs.  (16)  and  (17),  performing 
the  integration  over  (f>'  in  Eq.  (17),  and  comparing  terms  of  equal  order  in 
the  resulting  Fourier  series,  one  obtains  a  transfer  equation  for  each  Fourier 
component  of  the  Stokes  vector: 

dTm ( t  i/l 

u  a  \  ,U)  =1  2(t,u)  —  S™(t, u)  a  =  cots,  (29) 

dr 

where 

S?(t,«)  =  S£(t,U)  +  Q£(t,u)  (30) 

and  the  contributions  due  to  multiple  scattering  are 

s ?(T,u)  =  ^  j\u'\v^{T,u',u)l^{T,u'){l  +  80m) 

-Pr(r,u,,u)ir(T,u,)(l-<5om)}  (31) 

S T(r,u)  =  ^  j1  ^du'  +  50m) 

+  Pr(r,u,,u)I”l(r,u,)(l-<5om)}.  (32) 

While  the  cosine  modes  start  at  m  =  0,  the  sine  modes  start  at  m  =  1,  i.e. 
I ™(t,  u )  and  P™(r,  v! ,  u)  vanish  for  m  =  0. 

6.2.  Solution  of  the  RTE 

Half-range  quantities  in  slab  geometry 

In  slab  geometry  it  is  convenient  to  introduce  the  half-range  Stokes  vectors 

=  I+(r,0  <  tt/2,0)  =  I+(t,  (!,(/>) 

=  Ia(T,0  >  tt/2,</>)  =  (33) 
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where  a  =  c,  s,  and  fi  =  |w|  =  |  cos0|.  Note  that  the  sign  of  /z  is  given  by  the 
superscript  of  the  Stokes  vector.  For  the  Fourier  components  [see  Eq.  (22)],  we 
have  similarly 


I  ™+(T,d)  EE  I™+(T,<9<7r/2)=C+(T,/z) 

=  I%-(T,6>n/2)  =  I™-(T,n).  (34) 


6.2.1.  Formal  Solutions 

In  terms  of  these  half-range  quantities,  Eq.  (29)  may  written  as  (/i  =  |w|) 


dl™+{Fd) 

^  dr 
^  dr 


C+(r,/x)-S™+(r,M) 

ir(r,M)-s  r(r^). 


Using  an  integrating  factor,  we  readily  obtain  from  Eqs.  (35)  and  (36) 


(35) 

(36) 


-[C  +(r,M)e-TA1 


dl™+(r,/x)  lTTO+/_  ; 

dr  n  “  [T,N. 

Sq+(r,  /x)c-r/» 
d 

dI™~{T,n)  ,  lTm_,_ 

- -7 - k  -■!«  (U  d) 

dr  /i 


,~T/u 


cT/A* 


Cm— 

^ap 


(um)  t 


/M 


(37) 


(38) 


Homogeneous-slab  medium 

For  a  homogeneous  slab  with  constant  refractive  index,  we  adopt  the  following 
notation: 


1.  To  is  the  optical  depth  at  the  upper  boundary  (top  of  the  slab); 

2.  Tb  is  the  optical  depth  at  the  lower  boundary  (bottom  of  the  slab). 

Integrating  Eq.  (37)  from  to  r  and  Eq.  (38)  from  To  to  t  (t0  <  t  <  t^),  and 
solving  for  I”l±  (t,  /z)  we  obtain 


I f(^) 

=  I  ™+(Tfe,/x)e-^-^^  + 

[Tb 

(39) 

It  d 

i 

f-s 

(40) 

Jro  d 

Equations  (39)  and  (40)  show  that  for  cases  in  which  the  source  function  S™1*1  (t,  /z) 
is  known,  one  can  obtain  a  solution  to  the  radiative  transfer  problem  by  numer¬ 
ical  integration. 


Inhomogeneous  Slab  Multi-Layered  Medium 

The  vertical  variation  of  the  IOPs  in  a  slab  may  be  dealt  with  by  dividing  it 
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into  a  number  of  adjacent,  horizontal  layers  in  which  the  IOPs  are  taken  to 
be  constant  within  each  layer,  but  allowed  to  vary  from  layer  to  layer.  The 
number  of  layers  should  be  large  enough  to  resolve  the  vertical  variation  in  the 
IOPs.  In  such  a  multi-layered  medium,  consisting  of  a  total  of  L  layers,  we 
may  evaluate  the  integrals  in  Eqs.  (39)  and  (40)  by  integrating  layer  by  layer  as 
follows  (re -i  <  r  <  17  and  /i  >  0,  Tf,  =  tl,  to  =  0): 

[TL-S £+(t,/i)e-<*-T>/'‘  =  r  -S™+(f,/z)e^-^ 

Jt  A*  Jt  A*  ' 

+  E  I""  (41) 

n=t+ 

£ _ 

I"  -SZ~(t,n)e-[T~t)/tl  =  E  I  " 

Jo  M  ^  Jr„  .,  M 

/*r  dt 

+  /  -S£7(t,/i))e-<T-t>/'\  (42) 

Jr(_i  A* 

We  can  evaluate  the  integrals  in  Eqs.  (41)  and  (42)  either  numerically  or  an¬ 
alytically  if  the  source  function  is  known.  Analytical  integration  is  possible  in 
the  case  of  the  discrete  ordinate  method  to  be  discussed  below. 


6.2.2.  Discrete  ordinate  method 

Equations  (29)-(32)  can  be  rewritten  as: 

I ^  y'1idu'|pcm(r,u,,u)ir(r,u,)(l  +  ^om) 

-Pr(T,u»ir(T,u')}  -  Q?(t,u)  (43) 

I ?(r,u)  -  ^  j\u' 

+P ?(r,  u',  u )  I cm(r, «')}  -  Qsm (A  u)  (44) 

where  we  have  omitted  the  Kronecker-deltas  in  Eq.  (32)  and  the  second  Kronecker- 
delta  in  Eq.  (31)  since  the  cosine  modes  start  at  m  =  0,  whereas  the  sine  modes 
start  at  m  =  1,  i.e.  I™  and  P™  vanish  for  m  =  0.  The  discrete  ordinate  method 
consists  of  replacing  the  integration  over  u'  by  a  discrete  sum  by  introducing 
the  Gaussian  quadrature  points  Uj  (the  discrete  ordinates)  and  corresponding 


dl?(T,u) 

dr 


dr 
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weights  Wj.  One  obtains  for  each  Fourier  component: 


-Ed  £  "'i{(1+io”)P”<T'“''“‘)i"(T,“j) 

j  =  -N  ^ 

3  t^O 

(r,uJ-,ui)I™(T,'Uj)|  -  Q™(t,u;),  *  =  ±l,...,±iV  (45) 


-P 

(hui)  _  j mt  \ 

j  =  —  N  t 


+  P7*(r,uJ-,«i)I™(r)«J-)}  -  Qr(T.«i).  *=±1,...,±JV.  (46) 


The  convention  for  the  indices  of  the  quadrature  points  is  such  that  Uj  <  0  for 
j  <  0,  and  Uj  >  0  for  j  >  0.  These  points  are  distributed  symmetrically  about 
zero,  i.e.  U-j  =  — Uj.  The  corresponding  weights  are  equal,  i.e.  ui-j  =  Wj.  In 
a  plane-parallel  medium,  represented  by  a  horizontal  slab,  it  is  convenient  to 
consider  the  two  hemispheres  defined  by  u  <  0  (downward  radiation)  and  u  >  0 
(upward  radiation)  separately  and  to  make  us  of  the  half-range  quantity  /i  =  |rt| 
introduced  previously. 

It  should  be  emphasized  that  the  essence  of  the  discrete  ordinate  method  is 
to  convert  a  pair  of  coupled  integro-differential  equations  [Eqs.  (43)  and  (44)] 
into  a  system  of  coupled,  linear  differential  equations  [Eqs.  (45)  and  (46)]  by 
replacing  the  integrals  with  discrete  sums.  Combining  cosine  modes  of  the 
first  two  Stokes  parameters  (Jy  and  I±)  with  the  sine  modes  of  the  third  and 
fourth  Stokes  parameters  ( U  and  V)  in  Eqs.  (45)  and  (46),  one  obtains  a  set 
of  coupled  differential  equations  for  the  half-range  4IV  vectors  I™±(r, /x^)  = 
[T\\c  (T ,  ±/a),£C(t>  ±A*i)>  Vam(r,  ±m)]T,  i  =  l,...,N,  which  may  be 

written  in  matrix  form  as  (Weng,  1992:  Schulz  et  ah,  1999) 


-U =  (l  +  A?lc)  l?~  +  A-c  ir+  -  Qcm~ 

dTm+  .  -  , 

U-^y-  =  A-lcir  +  (1  +  A™C)I™+  -  Q™+ 

where  1  denotes  the  4 N  x  4 N  identity  matrix, 


Qcm± 


f  Q]fc(r-±Mi) 
OZc{T,±m) 
Qus(r,±^) 

V  Q^r.i/*) 


\ 

/  4iVxl 


(47) 

(48) 


(49) 
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and  U  is  the  4TV  x  4TV  diagonal  matrix 


(  \ 

U  =  u'i 

U, 

V  Uv 

where  the  TV  x  N  diagonal  submatrices  are  given  by 


(50) 


Uj_  =  Ujj  —  U[/  =  Uy  =  diag(/ui,  /^2,  ■  ■  • ,  Hn)T ■ 


The  matrices  A are  defined  in 

L  J  C. 

In  a  similar  manner,  by  combining  sine  modes  of  the  first  two  Stokes  pa¬ 
rameters  (I ||  and  /j_)  with  the  cosine  modes  of  the  third  and  fourth  Stokes 
parameters  ( U  and  V")  in  Eqs.  (45)  and  (46),  one  obtains  for  the  41V  vectors 

(r,  m)  =  K(r,  ±m):  I?S(T,  ±Mi),  U™(T,  ±m),  V™ (r,  ±^)}T ,  *  =  1 . JV, 

a  set  of  coupled  differential  equations,  which  may  be  written  in  matrix  form  as 
(Weng,  1992;  Schulz  et  al.,  1999) 


where 


-U 

u 


dl?- 

dr 


dr 


(l+A-jfr +Ar2sir+  -  Qlm- 
A™lsir  +  (i  +  A™j!r+  -  q;n+ 


Qlm± 


/  Q™s{t,±ih)  \ 

Q%c(Tt±m) 

V  / 


4iVXl 


(51) 

(52) 


(53) 


and  the  matrices  A"'-  are  defined  in 

lj  a 

Equations  (47)  and  (48)  or  Eqs.  (51)  and  (52)  may  be  re-written  as: 


or  as 


dVl 


dr 

dr 


\  m  j  m— 
Allot 

i  \m  jm+ 

'  A12a 

-  Qa' 

a  =  c,  s 

(54) 

\  m  jm— 
a21  a  La 

_L  A  m  Tm+ 

'  A22o:  La 

-  Q«  + 

Q 

II 

y* 

od 

(55) 

d 

t  m— 

/  a  m 
_  /  „  lla= 

Am  \ 
iyi2a  \ 

t  m— 

'  QS*" 

dr 

fm+ 

.  XQ! 

V  A m 

\  ^21a 

A  nn.  J 

^22a  / 

fm+ 

.  Aa 

.  Q™ +  . 

a  =  c,  s 


(56) 


where  A?"la  =  -IT^l  +  AftJ,  A?2a  =  -U"^,  A?la  =  U^A^, 
ASa  =  U-yi  +  A^2a),  and  Q£±  =  U^Q  am±. 

Combining  the  vectors  for  the  downward  and  upward  directions  by  defining 
the  8 TV  vectors  I™  =  and  Q™  =  [Q™-(W),  Q™+(Mi)F 
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(i  =  1, . . . ,  N)  we  may  write  the  coupled  system  of  differential  equations  in  the 
following  compact  form  (Weng,  1992;  Schulz  et  ah,  1999): 


(57) 


where  the  matrix  A™  is  given  by 


a  =  c,  s. 


(58) 


The  slab  is  divided  into  a  number  of  adjacent  layers,  large  enough  to  resolve 
vertical  changes  in  the  IOPs.  Equation  (57)  applies  in  each  layer  of  the  slab. 
As  described  in  some  detail  in  n  6.3,  the  solution  involves  the  following 

steps: 

1.  the  homogeneous  version  of  (57)  with  Q™  =  0  yields  a  linear  combination 
of  exponential  solutions  (with  unknown  coefficients)  obtained  by  solving 
an  algebraic  eigenvalue  problem; 

2.  analytic  particular  solutions  are  found  by  solving  a  system  of  linear  alge¬ 
braic  equations; 

3.  the  general  solution  is  obtained  by  adding  the  homogeneous  and  particular 
solutions; 

4.  the  solution  is  completed  by  imposing  boundary  conditions  at  the  top  and 
the  bottom  of  the  slab; 

5.  the  solutions  are  required  to  satisfy  continuity  conditions  across  layer  in¬ 
terfaces  in  the  slab; 

6.  the  application  of  boundary  and  layer  interface  conditions  leads  to  a  sys¬ 
tem  of  linear  algebraic  equations,  and  the  numerical  solution  of  this  system 
of  equations  yields  the  unknown  coefficients  in  the  homogenous  solutions. 

6.3.  Discrete  ordinate  equations  -  Compact  matrix  formulation 

As  shown  in  Section  6.2.2  each  Fourier  component  of  the  RTE  satisfies  the 
following  equations  (i  =  ±1, . . . ,  ±N): 


j  =  -N 
j^O 


(59) 


(60) 
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As  a  consequence  of  certain  symmetry  relations,  the  Fourier  components  of 
the  phase  matrix  have  the  following  form  and  van  dor  Mcc,  1983): 


/  pm 
/  -flic 

pm 

-*120 

0 

0  ^ 

pm 
-*  21c 

pm 

-*220 

0 

0 

0h_ 

II 

°  ^ 

0 

0 

pm 

-*33c 

pm 

-*340 

_  V  0 

PSJ 

V  0 

0 

pm 
-*  43c 

/  0 

0 

pm 
-*  13s 

Pus\ 

0 

0 

pm 
-*23  s 

pm 
-*  24s 

=  f  0 

pm 
-*  31s 

pm 
-*32  s 

0 

0 

11 

£3 

0  J 

Vpjl 

pm 
-*42  s 

0 

0  ^ 

(61) 


(62) 


where  we  have  introduced  (2  x  2)  block  matrices  for  notational  convenience. 


6.3.1.  “ Cosine ”  solutions 

As  discussed  in  Section  6.2.2,  by  combining  cosine  modes  of  the  first  two 
Stokes  parameters  (I ||  and  Jj_)  with  the  sine  modes  of  the  third  and  fourth 
Stokes  parameters  ( U  and  V)  in  Eqs.  (59)  and  (60),  one  obtains  a  set  of  coupled 
differential  equations  for  the  half-range  4 AT  vectors  I(”±(r,  fit),  which  may  be 
written  in  matrix  form  as 


dl?~ 

dr 

dl™+ 


dr 


\  m  t m—  \  \  m  ym+  _  rSm— 

^-11  cLc  '  a12 clc  Wc 

\  m  j m—  _i_  \m  ym+  _  (\rn-\- 

a21  cac  '  ^-22  cLc 


(63) 

(64) 


where  Ayjc  =  —  U 


-i/ 


-1  a  m  Am 
^-21 C5  ^-22c 


L/1  i  \m  \  A  m  _  tj-I  \m  \  m  _  tt- 
\±  '  Allc)i  Al$c  ~  U  A12c>  a21c  —  U 

U_1(l  +  A 22C),  and  Q™1*1  =  U_1Qcm±.  Here  U  is  a  4iV  x  4iV  diagonal  block 
matrix  [see  Eq.  (50)]  given  by 


U  =  diag(Uj_,  U || ,  Vv,  Uy)T 
where  the  N  x  N  diagonal  submatrices  are  given  by 

Uat  =  Uj_  =  XJ ||  =  U;/  =  XJV  =  diag(/xi,  /x2,  ■  ■  • ,  Hn)T 

and 


f  m± 


(  IFc(T,±IH)  \ 

/  0]pc(r,±Mi)  \ 

gm± 

QTc(r>  ±w) 

ujrir,  ±m) 

°c 

V  Kro(n±w)  y 

4iVxl 

V  Q^(t,±,ij;)  / 

(65) 

(66) 

(67) 


AN  x  1 


and  the  4 A/”  x  4 A  matrices  A™^  A^,  A^  and  A^  are  identical  expressions 
evaluated  at  different  pairs  of  angles: 


Lllc  V 


0  =  -' 


J  4 


Ai  +  <w)ps 

-prA 

l  ps 

ps  ; 

(68) 
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(69) 


A  72C(t) 

A  ?1c(t) 

ASe(r) 


Tu{t)  /(I  +  50m)Pllc 

j  4  (,  P2s 

ot(t)  /(I  +  <50m)Pic 
j  4  V  p2ms 

w(r)  /(l  +  <50m)P™ 
j  4  V  P2s 


prn 
rl  s 
■pm 
■*2c 


—  /Lti ,+Aij ;  z, .7=1,..., AT 


+/ii,—  //j;  i,j=l,...,N 


+fj,i,+Hj;  i,j=l,...,N 


(70) 


(71) 


where  we  have  used  u>_j  =  Wj. 

We  also  showed  in  Section  6.2.2  that  we  may  combine  the  vectors  for  the 
downward  and  upward  directions  in  Eqs.  (63)  and  (64)  by  defining  the  8 N  vec¬ 
tors,  i  =  =  [I™-(rt),I-+(rt)]T,  and  Q™  =  Q  ?+(^.)}T, 

so  that  we  may  write  the  coupled  system  of  differential  equations  for  these  “co¬ 
sine”  solutions  as 

dlm  ~  ~ 

'  =  A™I™  -  Q™  (72) 


where  the  matrix  A( 


dr 

is  given  by 


A  ™  = 


Lllc 
,  m 
l21  c 


L12c 


22c  /  SNx8N 


(73) 


We  arrived  at  Eq.  (72)  by  combining  cosine  modes  of  the  first  two  Stokes 
parameters  (I ||  and  I±)  with  the  sine  modes  of  the  third  and  fourth  Stokes 
parameters  ( U  and  V)  in  Eqs.  (59)  and  (60),  and  then  by  combining  the  resulting 
Eqs.  (63)  and  (64)  for  the  upper  and  lower  hemispheres  into  a  single  equation 
for  the  8 N  Stokes  vector. 

Remark:  if  the  incident  beam  has  no  U,  V  components,  as  for  the  disk- 
integrated  solar  pseudo-source,  then  the  homogeneous  “cosine”  solutions  ob¬ 
tained  by  solving  Eq.  (72)  completely  describe  the  homogeneous  solutions  be¬ 
cause  the  “sine”  solutions  described  below  vanish. 


6.3.2.  “Sine”  solutions 

In  a  similar  manner,  by  combining  the  sine  modes  of  the  first  two  Stokes 
parameters  (/ y  and  I±)  with  the  cosine  modes  of  the  third  and  fourth  Stokes 
parameters  (U  and  V)  in  Eqs.  (59)  and  (60),  one  obtains  for  the  4iV  vectors 
jm±  =  [J™(T)  ±m),  U™(t,  ikfii),  V™(t,  ±m)]T,  a  set  of  coupled  dif¬ 

ferential  equations  given  by  [see  Eqs.  (51)  and  (52)  and  Eqs.  (54)  and  (55)] 


drr 

dr 

dl™+ 


dr 


\  m  j m—  |  a  m  j m-\- 
-**■11  sLs  '  -**-12 sLs 

qt~ 

(74) 

\  m  j m—  a  m  jm+ 

-**■21  s  Ls  '  I*-22s 

-  Qr+ 

(75) 
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where  A™s  =  -IT^l  +  AftJ,  A 
U-^I  +  ASJ,  and  Q™±  =  U^Q 


m  T  j— 1  Am  \  m  t  7—  1  A  171  A  m  

\2 s  ~  u  ^12 s?  a21  s  —  U  **-21  ^22s  ~ 

m±  _  qqie  ma^rix  u  is  the  4j\r  x  4 A  diagonal 


matrix  given  by  Eq.  (65)  and 


(  ''>•*/'>)  \ 

/  Qfs(r,±W)  \ 

JT.(r,±w) 

Om±  — 

QTs(r,  ±lk) 

U™{t,  ±m) 

Vo is 

Quc(d  ±w) 

\  vcm(T,±m)  y 

4iVxl 

V  QVc(t,±im)  / 

(76) 


The  4iV  x  4iV  matrices  A")  s ,  Aj^s,  A™ls,  and  A^s  are  identical  expressions 
evaluated  at  different  pairs  of  angles: 


A  ?2a(r) 

A  ?ls(r) 


Mj) 

4 


(pm  Tym  \ 

“P£  (1  +  dom^TcJ 


(jym  Tym  \ 

^lc  rl,s  ] 

-P?s  (1  +  60m)P?c) 


(jym  Tym  \ 

*1  c  ^ls  ] 

-P%  (l  +  6om)P?J 


-IJ.i0.Hj-,  i,j=l,...,N 


+/j.j ,  —  fj,j ; 


(77) 


(78) 


(79) 


A  ?2s(t)  =  -Wj 


4  V- 


f  Pic 

Pii  A 

V~p£ 

(1  +  .WJP^y 

+fii,+Hj;  i,j=l,...,N 


(80) 


Combining  Ecis.  (74)  and  (75),  we  obtain  the  following  equation  for  the  8 N 


Stokes  vector 


dlm  -  ~ 

=  a^i  ?  -  Qs 


(81) 


which  is  of  a  form  identical  to  Eq.  (72).  Here 


A 


m 

s 


Am  \ 

**■123  \  .  t  m 

A  m  ’  $ 

A 22s /  8Nx8N 


Qs 


8JV  x  1 


The  “sine”  solutions  found  by  solving  Eq.  (81)  are  necessary  for  beam  sources 
with  non-zero  U,  V  components. 


6-4-  Discrete  ordinate  solutions 

The  vector  radiative  transfer  equations  [Eqs.  (72)  and  (81)]  have  the  same 
form  for  “cosine”  and  “sine”  combinations.  Thus,  we  may  rewrite  Eqs.  (72) 
and  (81)  as 

=  A ™(t)I-(t)  -  Q ™(t)  a  =  c,s.  (82) 
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6. 4-1.  Homogeneous  solution 

For  a  homogeneous  slab  with  a  constant  single-scattering  albedo  zu (r)  =  vo 
and  phase  matrix  P (t,u',u)  =  P (u',u),  the  matrix  A ™(r)  will  also  be  inde¬ 
pendent  of  position  r  in  the  slab.  Therefore,  we  may  seek  solutions  of  the 
homogeneous  version  of  Eq.  (82)  [setting  Q ™(r)  =  0]  of  the  following  form 

C(r)  =  G^e"^.  (83) 

Substituting  this  expression  into  Eq.  (82)  with  Q ”l(r)  =  0;  we  find 

A™G ™  -  -AmG™.  (84) 

Solving  this  algebraic  eigenvalue  problem,  we  get  8 N  eigenvalues  A"\ , A™8JV 
with  corresponding  eigenvectors  G”\, ...,  G™8iv.  The  general  solution  to  the 
homogeneous  version  of  Eq.  (82)  is  a  linear  combination  of  these  8 N  linearly 
independent  solutions  of  the  form  given  by  Eq.  (83). 


We  may  start  with  Eqs.  (63)  and  (64)  for  the  cosine  inodes  and  Eqs.  (74)  and 
(75)  for  the  sine  modes  and  rewrite  the  homogenous  version  of  these  equations 
as 


dr 


dr 

dl™+ 

dr 


\  m  t  m— 

Allot 

1  A m  fm+ 

'  A12a  La 

a  =  C,  S 

(85) 

\  m  j  m— 
a21  a  La 

1  a  m  fm+ 

'  ^-22 o: 

a  =  c,  s. 

(86) 

Dropping  the  superscript  m  and  the  subscript  a,  and  seeking  solutions  to 


Eqs.  (85)  and  (86)  of  the  form 

I ±(r, ±Mi)  =  G±e~kT  G±  =  G{k ,  ±m)  (87) 

one  finds  upon  invoking  the  symmetry  relations  A21  =  —  A12  and  A22  =  —An 
implied  by  the  symmetry  of  the  phase  matrix  (1  =  1, . . . ,  N ): 

— fcG~  =  An  G~  +  A12  G+  (88) 

-kG+  =  -A12G-  -  AnG+.  (89) 

Addition  and  subtraction  of  Eqs.  (88)  and  (89)  yields 

-k(  G-+G+)  =  (Au-Ai2)(G--G+)  (90) 

-fc(G"-G+)  =  (An  +  Ai2)(G“  +  G+)  (91) 


which  in  combination  gives 

(Aii-Ai2)(Aii+Ai2)(G-+G+)  =  k\ G"+G+)  (92) 

(Ah  +  Ai2)(Aii-Ai2)(G--G+)  =  k\ G--G+).  (93) 

Solving  the  eigenvalue  problem  given  by  Eq.  (92),  one  obtains  eigenvalues  k2 
and  corresponding  eigenvectors 

G+  +  G~  =  T  (94) 
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so  that  Eq.  (91)  yields 


G+-G-  =  }(An+A12)r. 

k 

Combining  Eqs.  (94)  and  (95),  one  finds 


G+ 

G~ 


r  +  —  (An  +  Ai2)r 
r  —  ^(A-n  +  Ai2)r 


which  on  noting  that 


G+ •  =  G+(—kj)  =  G~(kj)  =  G  J,  j  =  1,2,...,  AN 


(95) 


(96) 

(97) 


(98) 


so  that  G±(+kj)  =  GT(—kj)  provides  all  the  information  needed  for  the  solu¬ 
tion  of  Eqs.  (85)  and  (86). 

There  are  a  total  of  4 N  eigenvalues  =  k\,  i  =  1, . . . ,  4 N  and  hence  a  total 
of  8 N  eigenvalues  occurring  in  positive/negative  pairs: 

kj  =  v^>  0,  j  =  l,...,41V  (99) 

k-j  =  -kj,  j  =  1,...  ,41V.  (100) 


Since  both  real  and  complex  eigensolutions  may  be  present,  the  complex-variable 
eigensolver  DGEEV  of  LAPACK  [http://www.netlib.Org/lapack/lapack-3.5.0.html] 
may  be  used  to  get  the  solutions.  If  the  eigensolutions  are  known  to  be  real 
(as  for  Rayleigh  scattering),  the  faster  routine  ASYMTX  available  in  DISORT 

(Stamnes  et  ah,  1988,  2000)  may  be  used. 

The  solution  to  Eqs.  (85)  and  (86)  can  be  written  as  (*  =  1, . . . ,  N) 

4N  4  N 

I+(t)  =  I h(r,  +Hi)  =  E  CjGje-^  =  E[C,G+e-^  +  C-jGje+k‘T] 

J— — 41V  7  =  1 

3*  0 

(101) 

and 


4  AT  4N 

I-(t)  =  lh(r,  -/*)  =  E  CjGje-^  =  E  r/;,'  '  +  C.3G+e+kn 

J  — —  41V  7  =  1 

3*  0 

(102) 

or  more  compactly  as 


4iV 

Ifc  (r)  =  Mr,  ±Mi)  =  E^G fe~k>T  +  CL,- Gje+M  (103) 

i=i 


where  { Cj }  and  {C-j}  are  constants  of  integration  to  be  determined  by  the 
boundary  conditions. 
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Since  both  the  eigenvalues  {kj}  and  the  eigenvectors  G±{kj)  =  GT(- kj) 
may  be  complex,  it  is  desirable  to  rewrite  Eqs.  (101)  and  (102)  in  terms  of 
real  quantities.  To  that  end,  one  may  distinguish  between  real  and  complex 
eigenvalues  {kj}.  The  eigenvalues  (and  the  associated  eigenvectors)  of  a  real 
matrix  are  either  real  or  they  occur  in  complex  conjugate  pairs  Strang,  2005). 
Thus,  letting  Jr  denote  the  number  of  real  eigenvalues  and  Jx  denote  the  number 
of  complex  conjugate  pairs  of  eigenvalues,  one  may  write  Eqs.  (101)  and  (102) 
as 

=  +  (104) 

where  the  real-valued  solutions  are 

Jr 

=  £[QG+e-^T  +  C-jG~  e+kjT]  (105) 

i=i 

Jr 

KJr)  =  ERG-e-^  +  C-j  G+e+kn  (106) 

j= i 

and  the  complex-valued  solutions  are 


J , 


=  Vr,,,G;v  +  cx, 

3  =  1 

(107) 

-c* 

1  H 

hH 

=  £,[C*,jGZije-k*^  +  Cx,-3  Gttje+k-n 

7=1 

(108) 

where  the  subscript  x  (last  letter  in  the  word  complex)  is  used  to  indicate  that  a 
quantity  is  a  complex  number  or  eigenvector.  Taking  the  real  parts  of  Eqs.  (107) 
and  (108),  we  find  {Cxd  =  CRd  -) -iCij) 

mtnir)} 

Jx 

=  E  {6^F+(r,  kj)  +  Cr^Fr{t,  -kj) 

j= i 

~  (r,  kj)  —  Ci-jFj  (r,  — fcj)  j 

(109) 

=  E  {^fldF/f.(r:  kj)  +  CR_jFR(r,  —kj) 

3=1 

-  CIdFj{r,kj)  -  G^jF+i^-kj)} 

(HO) 

where 

F  %{r,±kj) 

=  ;^G;.?j^{(  •-)  ^{G;Jd{f  ;A'-| 

(111) 

F  f(r,±kj)  = 

[5R{G±  .}9{e=F*-.J-}  +  3{G ^{eTl-iT}]  • 

(112) 

Note  that  all  the  8 N  =  2{Jr  +  2JX)  constants  (2 Jr  values  of  C±j,  and  2 Jx  values 
each  of  CR.±j  and  Cit±j)  in  Eqs.  (105)-(112)  are  to  be  determined  by  applying 
the  boundary  conditions. 
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6. 4-2.  Vertically  inhomogeneous  medium 

For  a  vertically  inhomogeneous  medium,  one  may  divide  the  slab  into  L 
adjacent  homogeneous  layers  labeled  by  the  index  i  =  1  The  IOPs 

(the  phase  matrix  and  the  single-scattering  albedo)  are  assumed  to  be  constant 
within  each  layer,  but  they  are  allowed  to  vary  from  layer  to  layer  as  required  to 
adequately  resolve  their  vertical  variation.  Then,  the  solution  arrived  at  above 
[Eq.  (83)]  will  be  applicable  to  each  layer  separately.  The  exposition  may  be 
simplified  by  rewriting  the  homogeneous  solution  given  by  Eq.  (103)  for  mode 
a  =  c  or  s  and  layer  denoted  by  £  as 

AN 

Y,\-C-^G-%ee+k^lT  +  QG^e'Cj,T]-  (H3) 

j= i 

Defining  !^(r)  =  [I"7(r),  I™+(r)]T,  Caje  =  [C-ajt,  Caje]T,  Gaje  =  [G-aje,Gaj 
one  may  write  Eq.  (113)  representing  the  solution  to  the  homogeneous  part  of 
the  RTE  for  a  Fourier  mode  denoted  by  the  subscript  a  and  a  component  de¬ 
noted  by  the  superscript  m  in  a  layer  denoted  by  subscript  l  more  compactly 
as 

8N 

C/V)  =  E  C^yG^e-*^  £  =  1, . . . ,  L.  (114) 

3= 1 


a£  (T 


6. 4-3.  Particular  solution 

The  source  term  Q™(r)  in  Eq.  (82)  for  the  cosine  modes  is  given  by 


n-  =  ( Q“” 

-  l  Qm+ 


Qcm±  = 


8JVxl 


f  Qfc(r,±W)  ^ 
QZc{t,±ih) 
QV,(T,±tH) 

V  Q$b(t,±ih)  J 


4Nx  1 


and  for  the  sine  modes 


Am  =  ( QT~ 

Q s  ~  (qt+ 


Qln±  = 


8N  x  1 


(  Qfs(T,±rt)  \ 
QZ(p  ±m) 
QVc(r,±m) 

V  Qyc(T,±m)  ) 


4Nx  1 


where  ( a  =  c,  s)  Q  Q±('r,  Mi)  =  U  1QcJ?1±(t,  /r*)  and 

Q«m±fo  w)  =  X™±(r,/xl)  e~T^°  +  5ca  B0m(r). 


(115) 


Here  X”l±(r,  m)  =  P ™(r,  -/x0,  i/q)  S6.  Thus,  the  source  term  Q ™(r,  /Zj) 

in  Eq.  (82)  can  be  written  for  each  layer  denoted  by  subscript  £\ 

Q S(Mi)  =  [XS(w)  e~T^°  +  5ca  B0m,,] 


(116) 
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where  X£(w)  =  U^X™-  (m),  X™+(Mi)]T,  and  B0m>*  =  <50ro  [1  -  w/]  U"^ 
[see  Eq.  (28)]. 


Cosine  modes 

Combining  the  two  upper  rows  of  the  matrix  P™  and  the  two  lower  rows  of  the 
matrix  P™  given  by  [see  Eqs.  (61)  and  (62)] 


/  pm 
/  -'lie 

pm 
-*12  c 

0 

0  ) 

/  o 

0 

pm 

-*13s 

pm  \ 
-*14  s  \ 

pm 
-*21  c 

pm 
-*22  c 

0 

0 

.  -pm  _ 

0 

0 

pm 

-*23s 

pm 
-*24  s 

0 

0 

pm 
-*33  c 

pm 

■*34c 

5  -*-3 

pm 

-*31  s 

pm 
-*32  s 

0 

0 

V  0 

0 

pm 
-*43  c 

pm 
-*44  c) 

,  pm 
\-*41s 

pm 
-*42  s 

0 

o  / 

one  obtains  a  matrix  P™  and  a  corresponding  vector  S;,  defined  as 


(Pile 

pm 
-*12  c 

0 

pm 

-*21c 

pm 

r22c 

0 

0 

pm 
-*31  s 

pm 
r32  s 

0 

0 

,  pm 
Y*41s 

pm 
-*42  s 

0 

<V 

8  iVi  x8iVi 

(S  ||6  \ 

s±b 

Sub 

\SvbJ  8ArlXl 


(117) 


Hence,  for  the  cosine  modes  we  obtain 


X 


m 

c 


tz7(r) 

/  pm 
'  -'lie 

pm 
-*12  c 

0 

pm 
-*21  c 

pm 

r22c 

0 

0  I 

47T 

pm 
■*3  Is 

pm 

-'32s 

0 

0 

Sub 

pm 
-*42  s 

0 

8JVi  x8ATi 

\SvbJ 

™(T) 

47 r 


prs  h 


p m 

4t r  c’h 


(118) 


where  the  vector  P™b  =  P™Sb  is  given  by 


pm 

*6,  c 


/P(?c(r,  -Mo,  ~M?)%  +  "Mo,  -M?)£lA 

P%c(t,  -Mo,  -M?)S||6  +  ^c(r,  -Mo,  ~tf)S±b 
PR,(t,  -Mo,  ~M“)%  +  P^(t,  -Mo,  -M?)£l6 
^41S(P  -Mo,  — Mf )'S'||6  +  P42S(p  -Mo,  ~Hi)S±b 
Pi ic(r,  -Mo,  +m“)%  +  ^I2c(r,  -Mo,  +M?)5’_Lb 
A7 ic(r,  -Mo,  +M?)%  +  AT2c(r>  -Mo,  +M?)>Su, 
^3is(r,  “Mo,  +M?)%  +  P^a(r,  -Mo,  +M?)'S'_Lb 
\p4is(r,  “Mo,  +Mi  )^||fc  +  P42S(r,  “Mo,  +Mf)'S,J-fe/ 


8Ni  x  1 


Sine  modes 

For  the  sine  modes  one  may  proceed  in  a  similar  manner  by  combining  the  two 
upper  rows  of  the  matrix  P”1  and  the  two  lower  rows  of  the  matrix  P™  to  obtain 


(0 

0 

pm 
■*1  3s 

pm  \ 
-'14s  \ 

0 

0 

pm 
-*23  s 

pm 
-*24  s 

0 

0 

pm 
-*33  c 

pm 
-*34  c 

\o 

0 

pm 
-*43  c 

pm  , 
-*44  c/ 
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Hence,  for  the  sine  modes,  one  obtains 


p.r 

s  _  47T  s 


s,  =  5Mp- 


47T 


where  the  vector  P™fc  =  P( 


Sb  is  given  by 


pm  _ 

s,b 


/P^s(T,-M0,~K)SUb 

P2™s(T,-M0,-K)SUb 

p%c(T,-Mo,-K)sUb 

P43c(T’-V0,-K)SUb 
Pl3s(T,~/M),+Vi)Sub 
Mo,  +  Mi)  Sub  ■ 
Mo,  +Mi)Sub  - 
(t,  ~Mo,  +Mi)SUb  - 


13s  V 
Dm  / 
23s  V 


33  c(T, 


P, 

V^43c 


Pns(r,-Mo,-Mt)SVb\ 

P24S{r,-Ho,-Mi)Svb 
^34c(r)  —MO,  ~Mi)Svb 
PaIc{t,  ~Mo,  —Mi)Svb 
■  Pusir,  -mo,  +Mai)Svb 
P2ls(T,~Mo,+Mi)Svb 
PSlc(T,-Mo,+Mt)Svb 
P4Ac(T,-M0,+Mi)SVbJ 


(120) 


8N±  x  1 


From  this  expression  it  follows  that  if  the  incident  beam  source  has  no  U  and 
V  components,  then  the  sine  modes  of  the  (I\\,I±)  and  the  cosine  modes  of  the 
( U ,  V)  Stokes  vector  components  vanish,  in  which  case  Eq.  (72)  provides  the 
complete  homogeneous  solution  as  I™  =  0. 


Beam  source 

Consider  now  solutions  for  each  of  the  source  terms  in  Eq.  (116).  Thus,  in  a 
layer  denoted  by  subscript  i,  for  the  term  e~T^°,  where  is  assumed  to 
be  constant  within  layer  £,  consider  a  particular  solution  of  the  form 

!sw  =  ZSe-^o.  (121) 

Substitution  in  Eq.  (82)  with  Q™  =  e~T^°  leads  to  the  following  system 

of  linear  algebraic  equations  [1  is  the  8N  x  8N  identity  matrix] 

[A™  +  (-)l]Z^  =  XS  (122) 

Mo 

which  can  be  solved  to  yield  a  particular  solution  vector  Z™(. 

Thermal  source 

From  the  thermal  source  term,  there  is  a  contribution  only  for  cosine  mode  m  = 
0  in  the  Fourier  expansion;  there  is  no  contribution  for  any  of  the  sine  modes.  To 
get  an  approximate  particular  solution  to  the  inhomogeneous  radiative  transfer 
equation,  one  may  make  polynomial  approximations  to  the  Planck  function  and 
the  source  vector  S t(r): 


I< 


B(T(r))  =  J2b^k 

k= 0 

K 

Qf  =  U-!(1- W(r))^Dfcri 


k=0 
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where 


Di,  = 


/Dlfc\ 

Dofc 

Difc 

VDofcy 


Difc  = 


Dofc  = 


(123) 


8iVx  1 


2N  x  1 


2Nx  1 


where  it  has  been  taken  into  account  that  the  thermal  source  is  unpolarized. 
Since  the  thermal  source  term  is  non-vanishing  only  for  the  cosine  mode  m  =  0, 
Q™  =  SomQt,  Eq.  (82)  becomes: 


di?(r) 

dr 


=  A™(t)I™(t)  -  S0mQt(T) 


Thus,  by  substitution  in  Eq.  (124)  the  assumption: 


K 


I?(r)  =  50mY,X 


t  k 


k—0 


one  obtains 


K 


K 


K 


S0m  J2  k  XfeTfe_1  =  60mA™  J2  -  5om(  1  -  w(r))  ^  D krk 

k= 1  k—0 

r  K+l  K+l 

=  ^Om 


k—0 

a™  53  n^-1  -  (i  -  w{t))  53 

*-  k= 1  fe=l 

Comparison  of  terms  having  equal  powers  of  r  leads  to 

k  <  K  +1  :  k  X*  =  A°X*  -  (1  -  tn(r))Dfc_1 

/c  =  A  +  l  :  0  =  (A°X^-  (1-  w(r))Dx. 

The  first  order  approximation  (A'  =  1)  is  usually  sufficient  and  yields 

and  then  we  find  Xq  from 

A°X*  =  (1  -  tu(t))D0  +  X*. 


(124) 


(125) 


(126) 

(127) 


6. 4-4-  General  solution 

For  each  Fourier  component  denoted  by  the  superscript  m,  the  general  solu¬ 
tion  is  a  combination  of  the  homogeneous  solutions  and  the  particular  solutions 
for  beam  and  thermal  sources.  Thus,  the  general  solution  is 

8  N 

W  =  53C™.,G^e-fe^  +  ZSe-^ 
i= i 

+  <5om<5ac[Xg ^  +  X*  ^r],  £=1,...,L.  (128) 
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Since  each  of  the  solution  vectors  in  the  upper  slab  are  of  dimension  8 N,  for 

example  G%t  =  fe(l),  G%t(2), . . . ,  G™e(8N)]T  =  =  1 . 8 N, 

where  the  lower  case  letter  g  has  been  used  to  denote  the  vector  component, 
one  may  rewrite  Eq.  (128)  in  component  form  as  follows  (i  =  1, . . . ,  8 N) 

8N 

E  ^&i)e-k^T  +  CA*)e~T/M 

j= i 

^Om^ac (^)  ^  —  I?  •  •  •  >  A.  (129) 

6.Jf.5.  Boundary  conditions 

To  complete  the  solution  one  must  apply  boundary  conditions  in  order  to 
determine  the  unknown  coefficients,  the  C™/ s,  in  Eq.  (128).  Mathematically, 
one  is  faced  with  a  two-point  boundary  value  problem,  which  requires  a  spec¬ 
ification  of  the  radiation  field  incident  at  the  top  and  the  bottom  of  the  slab. 
In  addition,  at  layer  interfaces  within  the  slab,  the  radiation  field  is  required  to 
be  continuous  because  the  refractive  index  is  assumed  to  be  constant  within  the 
slab.  So,  in  summary,  the  following  conditions  are  required 

1.  At  the  top  of  the  slab  the  incident  Stokes  vector  must  be  specified; 

2.  The  Stokes  vector  must  be  continuous  across  layer  interfaces  within  the 
slab,  I at+iin)  =  iaein); 

3.  At  the  bottom  boundary  the  Stokes  vector  must  be  specified  in  terms  of 
the  bidirectional  polarized  (reflectance)  distribution  function. 

Implementation  of  these  conditions  leads  to  a  system  of  linear  equations,  and  its 
solution  yields  the  desired  coefficients,  the  C'™/s.  For  each  layer  there  are  8N 
equations  corresponding  to  the  Stokes  vector  components  at  the  8N  quadrature 
points.  For  L  layers  one  obtains  8 N  x  L  equations.  Hence,  one  obtains  a  linear 
system  of  equations  of  dimension  (8 IV  xL)x  (8 N  x  L) ,  and  the  solution  of  this 
linear  system  of  equations  will  determine  the  (8 N  x  L)  unknown  coefficients, 
the  C™/s,  required  to  complete  the  solution. 
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